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ABSTRACT: 

A quadrilateral, eight-node, mixed plate element based upon a recent 3D, variable kinematics zig-zag 
plate model by the authors is developed. Representation can be refined across the thickness though 
the number of unknown does not depend on the number of constituent layers, the nodal variables being 
fixed. The out-of-plane stress continuity is a priori satisfied at the interfaces and the material 
properties are allowed to step-vary moving over the plane of the structure. In order to obtain a CO 
model, all derivatives of unknowns are converted with a technique that does not introduce additional 
d.o.f. Equilibrium equations are satisfied just in approximate integral form, mid-plane displacements 
and interlaminar stresses, the nodal d.o.f, being interpolated using the same standard serendipity 
shape functions. As shown by numerical results, this does not result in any precision loss. Accuracy, 
solvability and convergence are assessed considering multi-layered monolithic and sandwich-like 
structures with different boundary conditions and abruptly changing material properties across the 
thickness, for which exact solutions are available. Also a bonded joint is studied, which is treated as a 
laminate with spatially variable properties. In all these cases, the element is shown to be fast 
convergent and free from spurious, oscillating results. 

KEYWORDS: partial hybrid formulation; intra-element equilibrium in an approximate integral 
form; fixed d.o.f; hierarchic representation; CO plate element; 



I. INTRODUCTION 

Laminated and sandwich composites increasingly find use, as they enable design and construction of 
structures that achieve the target requirements with a lower mass than their metallic monolithic counterparts and 
they also offer the possibility to optimize the structural performances by properly choosing the fibre orientation 
and the stacking lay-up (see, e.g. Sliseris and Rocens [1]). 

Manufacturing technologies such as automated fibre -placement (Barth [2]) pave the way to the spread 
of variable-stiffness composites in which the fibres follow curvilinear paths. Variable strength, stiffness and 
dissipation properties offer many advantages, as shown e.g. in Refs. [3] and [4], since the constituent materials 
can be designed using advanced optimization techniques for keeping maximal the overall performance 
parameters. At the same time, it is possible to relax the critical local stress concentrations inherent to the 
inhomogeneous microstructure of composites that give rise to damage. A recent, comprehensive discussion of 
the mechanisms of damage formation and evolution of composites and of their modelling is given by Cardenas 
etal. [5]. 

Since the local damage can have very harmful effects, the simulations should accurately account for the 
stress fields responsible of the damage rise and growth, their strong concentrations at the interfaces, their effects 
on strength, stiffness and durability and the warping, shearing and straining deformations of the normal 
produced by the material discontinuities and the much bigger elastic moduli and strengths in the in-plane 
direction compared to those in the thickness. These so called zig-zag and layerwise effects should be described 
into details but with an affordable computational effort, as discussed by Chakrabarti et al. [6], Qatu et al. [7] and 
Zhang and Yang [8]. 

Displacement-based models and related finite elements should feature appropriate discontinuous 
derivatives of displacements in the thickness direction at the interfaces, while mixed models should assume 
appropriate continuous, self-equilibrated stresses. Composite plate and shell theories and elements satisfying 
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these requirements have been developed either in displacement -based or mixed form using different approaches. 
Displacement-based finite elements require fine meshes and a large computational effort to accurately determine 
stresses, since they are obtained from displacements (constitutive equations) or from in-plane stresses 
(integration of local equilibrium equations), thus accuracy deteriorates. Mixed/hybrid models and elements offer 
the advantage of a more easily enforcement of interfacial and boundary constraints, the stresses being treated as 
primary variables separately from displacements. Moreover, they provide more accurate results than 
displacement-based models with the same meshing and with a comparable processing time, as shown in the 
literature. On the other hand stability, convergence and solvability are more complex than those of their 
displacement-based counterparts. 

An extensive discussion of the various techniques used to account for the layerwise effects and 
extensive assessments of their structural performances are presented by Chakrabarti et al. [6], Matsunaga 
[9]Error! Reference source not found., Chen and Wu [10], Kreja [11], Tahani [12] and Gherlone [13]. 
Accuracy of finite element models for layered and sandwich composites is assessed by Carrera et al. [14]. Finite 
element models are also discussed in the papers by Chakrabarti et al. [6], Zhang and Yang [8], Shimpi and 
Ainapure [15], Elmalich and Rabinovitch [16], Dau et al. [17] and for what concerns mixed/hybrid methods in 
the book by Hoa and Feng [18], and in the papers by Feng and Hoa [19], Desai et al. [20], Ramtekkar et al. [21] 
and To and Liu [22]. Recent examples of finite elements for analysis of composites are those by Zhen et al [23], 
Cao et al. [24] (hybrid formulation) and Dey at al. [25]. 

Accuracy and computational costs of layerwise models and elements can considerably vary, depending 
on their number of unknowns. A separate representation in any computational layer results in a high accuracy 
(see, e.g. Reddy [26]), but since the number of variables increases with the number of physical/computational 
layers, computational costs can become unaffordable when analysing structures of industrial complexity. 
Models based on a combination of global higher-order terms and local layerwise functions are often used for 
analysis of composites, as they have been proven to be equally accurate with a much lower computational effort 
(see, e.g. Elmalich and Rabinovitch [16] and the references therein cited). The zig-zag models inspired by Di 
Sciuva [27] have also found many applications and related finite elements are widespread, being low cost partial 
layerwise models with just the three mid-plane displacements and the two transverse shear rotations as 
functional d.o.f., like equivalent single-layer models. Physically-based continuity functions are incorporated in 
the displacement field in order to a priori satisfy the continuity of out-of-plane stresses at the layer interfaces, 
whose expressions are determined once for all. Initially, just the piecewise variation of in-plane displacements 
across the thickness was considered, the transverse displacement being assumed constant. Refinements have 
been brought over years in order to incorporate a variable transverse displacement, since it has a significant 
bearing for keeping equilibrium in certain cases. A great amount of research work was also done to successfully 
treat panels with low length-to thickness ratio and abruptly changing material properties, like sandwiches, which 
can be described as multi-layered structures whenever a detailed description of local phenomena in the cellular 
structure is unnecessary [28]. Sublaminate models having top and bottom face d.o.f. were developed in order to 
stack computational layers [29]. The displacement field was recast in a global-local form to accurately predict 
stresses from constitutive equations (see, Li and Liu [30], Zhen and Wanji [31]), though at the expense of a 
larger number of functional d.o.f., because post-processing operations are unwise for finite elements and cannot 
always be precise [32] . 

As a result of these refinement studies, to date zig-zag models and elements offer a good accuracy with 
the lowest computational burden. Regrettably, derivatives of the functional d.o.f. are involved by such 
physically-based zig-zag models, and consequently they should appear as nodal d.o.f. in the finite elements, 
implying to use CI or high-order interpolation functions with a lower computational efficiency than COones. 
Techniques such those proposed by Zhen and Wanji [31] and Sahoo and Singh [33] could be employed for 
converting derivatives, but they result in an increase of the nodal d.o.f and thus of the memory storage 
dimension. The Murakami's zig-zag function just based upon kinematic assumptions is often used as local 
layerwise function since it overcomes the problem and it is much easier to implement than the physically 
consistent zig-zag functions, though it is not always equally accurate. The assessments carried out by Gherlone 
[13] proven that Murakami's zig-zag is accurate for periodical stack-ups, but not for laminates with arbitrary 
stacking sequences, or for asymmetrical sandwiches with high face-to-core stiffness ratios, like when a weak 
layer is placed on the top or bottom to simulate a damaged face. 

Aimed at carrying out the analysis of multi-layered and sandwich composites having abruptly 
changing material properties with the minimal computational burden, the recently developed physically-based 
zig-zag model of Ref. [34] was developed assuming a through-the-thickness variable representation of 
displacements that can be locally refined like for discrete layer models, though its functional d.o.f. are fixed ( 
the classical displacements and shear rotations of the normal at the mid-plane). A high-order piecewise zig-zag 
variation of all the displacements was assumed in order to account for the core's crushing behaviour of 
sandwiches, for keeping equilibrium at cut-outs, free edges, nearby material/geometric discontinuities and to 
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predict stresses caused by temperature gradients, because in these cases the piecewise variation of the 
transverse displacement is as important as that of its in-plane counterparts (see, e.g. [23] and [35]). The physical 
constraints, i.e. the interfacial stress contact conditions, stress boundary conditions at the bounding faces and the 
indefinite equilibrium equations at selected points across the thickness were satisfied without increasing the 
number of unknown variables because symbolic calculus was used to obtain automatically and once for all 
relations in closed-form as functions of the assumed five d.o.f. Accurate results were obtained in numerical 
applications with a computational effort considerably lower than for other layerwise models, thus model [34] is 
of practical interest toward development of finite elements. 

In the present paper, the energy updating technique of Refs. [36] - [38], here referred as SEUPT, is 
used to convert derivatives without resulting in an increased number of nodal variables. SEUPT was originally 
developed as an iterative post-processing technique for improving the predictive capability of shear deformable 
commercial finite plate elements. The present revised version of SEUPT consists of a technique that allows 
obtainment of an equivalent CO version (EM model) of the zig-zag model [34] by the energy standpoint, which 
is used to develop an efficient plate element in mixed form having the displacements and the interlaminar 
stresses at the reference mid-plane as nodal d.o.f. A priori corrections of displacements are obtained in closed- 
form through SEUPT using symbolic calculus by making the energy of the model [34] free from derivatives 
equal to that of its counterpart containing all the derivatives. Hermite's polynomials are assumed in the energy 
balance equations for the consistent zig-zag model, while Lagrange's polynomials are used for the equivalent 
counterpart model free from d.o.f. derivatives. The equivalent CO displacement fields and the out-of-plane 
stresses derived in this way are used to develop the present mixed element. It is a partial hybrid one, because, 
according to Ref. [19], it only includes out-of-plane stresses as nodal d.o.f., just these stresses having to be 
continuous across the material interfaces and being essential for analysis of composites. The stresses by the EM 
model are modified adding new continuity functions that restore the equilibria at the interfaces, because stress 
continuity could be lost in the equivalent model. 

Displacement and stress interpolating functions all of the same order are considered, the same CO 
interpolation by standard serendipity polynomials being used for all the nodal d.o.f. , so the intra element 
equilibrium equations are satisfied just in an approximate integral form. This option was chosen because it 
reduces the effort for developing the element without losing accuracy, as shown in the pioneering paper by 
Loubignac et al. [39] and in the book by Nakazawa [40]. 

Accuracy, solvability and convergence behaviour of the present mixed element will be assessed 
comparing its results for laminates and sandwiches with different boundary conditions and abruptly changing 
material properties with the exact solutions computed with the technique by Pagano [41]. The sample cases by 
Zhen et al. [23], Vel and Batra [42], Icardi [43], Brischetto et al. [44], Tessler et al. [45] and Robbins and Reddy 
[46] are considered. As an application to a practical case, the stress field of a double lap joint computed by the 
present mixed element is compared to the exact solution by Nemes and Lachaud [47]. 

The paper is structured as follows. First, the basic features of the zig-zag plate model are overviewed. 
Then, SEUPT is recast in a form suited for the present element. Finally, the element is developed and applied to 
the analysis of sample cases. 



A rectangular Cartesian reference frame (x, y, z) with (x, y) on the middle surface of the plate (£2) and z 
normal to it is assumed as reference system. The displacement field is assumed as the as the sum of four 
separated contributions: 

u(x, y, z) = U ° (x, y, z) + U ' (x, y, z) + U ° (x, y, z) + U °' v (x, y, z) 

v (x , y , z) = V ° (x , y , z) + V ' (x , y , z) + V ' (x , y , z) + V ' ~ ''' (x , y , z) ^) 
w(x, y , z) = W " [x, y , z) + W ' [x, y , z) + W c (x, y , z) + W (x,y,z) 

Here the symbols u, v and w respectively define the elastic displacements in the directions x, y and z. 

The first three contributions are the same of the model developed in Ref. [34], the latter one was added 
in Ref. [38] to treat properties and/or lay-up suddenly varying in x, y directions, like across adherends and 
overlap of bonded joints. 

2.1 Basic contribution □ 

The basic contributions, indicated with superscript 0 or in compact form as A 0 , repeat the kinematics of 
the FSDPT model, thus they contain just a linear expansion in z: 



II. 



STRUCTURAL MODEL 




(2) 
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W (x,y,z) = w (x,y) 

Eq. (2) is well known and it has been used in many applications, therefore it does not need further 
explanations. Here it is only reminded that the displacements u°(x, y), v°(x, y), w°(x, y) and the transverse shear 
rotations at the middle plane y x °(x, y) and y y °(x, y) represent the five functional d.o.f. 

2.2 Variable kinematics contribution □' 

The contributions with superscript ', indicated in compact form with the symbol A 1 , allow the model 
adapting to local variation of solutions, enabling a refinement across the thickness. Therefore these contributions 
have a representation that can vary from point to point across the thickness: 

U'(x,y,z) = A xl z+ A^z 2 + A t3 z 3 + A i4 z" + ... + A a z" 
V ' (x, y, z) = A yl z + A y2 z 2 + A y3 z 3 + A v4 z* + ... + A yn z" 

W\x,y,z) = A zl z + A :2 z 2 + A_ 3 z 3 + A z4 z* + ... + A :n z" 

The unknown coefficients A xl . . . A zn are computed by enforcing the conditions: 

a I" = 0 a I, = 0 

(4) 

* ° " ' ° (5) 

a 1 = p 1 °" Jl = p '/ ( 6 ) 

a I" = 0 a 1=0 

(7) 

and the equilibrium at discrete points across the thickness: 

a a t + a x + a xz z =0 

a + a + a =0 

xy ,x yy>y yz,z 

a + a + a =0 

(8) 

The symbols (k) z + and (k) z" indicate the position of the upper + and lower " surfaces of the kth layer, while 
the superscript(k) designates the quantities that belong to a generic layer k. As customarily, a comma is used to 
indicate differentiation. Please notice that the expressions of the unknowns A xl ... A zn are computed through 
automatic, symbolic calculus (MATLAB® symbolic software package) without any increase in the number of 
unknowns, being obtained in closed-form as functions of the d.o.f. and of their derivatives. Contributions (3) 
neither results into a considerably larger computational effort, nor into a larger memory storage dimension, 
because irrespectively of the number of points chosen, the number of d.o.f. is fixed because the constraint 
relations are obtained in closed-form once for all just in terms of the variables u°(x, y), v°(x, y), w°(x, y), y x °(x, 
y), y y °(x, y) and of their derivatives. Because derivatives are unwise for the development of finite elements, they 
will be converted with the technique described forward. 

2.3 Zig-zag piecewise contribution A c 

The piecewise terms Ac, which constitutes the zig-zag contributions, are assumed as follows: 

U °{x,y) = 2 <K*(*,yXz - z k )H t + £ C* (x,y)H t 

v ' (x, y) = x * ' (*• yx z - i t ) H t + Z c ' (*> y~> H t 

(9) 

W "(x,y) = ^f " (x,y)(z- z t )H t + £ £2 * (x, y )( z - z k )* H k +•£ C*(x,y)H k 

The terms of Eq. (9) make the displacements continuous and with appropriate discontinuous derivatives 
in the thickness direction at the interfaces of physical/computational layers, thus allowing to a priori fulfil the 
continuity of interlaminar stresses at the material interfaces. In details, terms <D x k , <D y k are incorporated in order 
to satisfy: 



while *T k , O k are computed by enforcing the following continuity conditions: 



(10) 
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Z Z 

,* " r (11) 
which directly derive from the local equilibrium equations as a consequence of the continuity of 
transverse shear stresses. 

The continuity of displacements at the points across the thickness where the representation is varied is 
restored by the functions C u k , C v and C w \ whose expressions are computed imposing: 



(12) 



2.4 Variable in-plane representation A CJP 

Differently to [34]Error! Reference source not found., it is supposed that the in-plane representation 
of displacements can vary, consequently to a sudden change of the material properties moving along the x or y 
direction. To this purpose, contributions Ac ip are incorporated in a form similar to that used in [38]: 

ST ST 

W = II X'UjKx-^JS. + n y y (x,y)(y-y k )H k + 

7 = 1 k=\ ./ = 1 k=\ 

ST ST 



/=1 k = 1 
S T 



(13) 



"'■'=11 XU.jK'-^lfl. + n le;{x,y)(y-y k )H t + 

7=1 k=l 7=1 k=\ 

ST ST 

XX Xu.yX^-VX+ZE X{x,y){y-y k ?H k+ ... 

j = \ k = \ j = \ k=\ (13)' 

In Eqs. (13) and (13)' the exponent of (x - x k ) n , (y - y k )" is chosen in order to make continuous the 
stress gradient of order n. For instance, the continuity conditions a xx \(k) X + = c.o-W-, ffyyWv+ = c^W- and <7 a J(jy I+ = 
°zzW- require just first order terms of (x - x k ), while the continuity of gradients a XXiX \(k) X + = OxxJmx-, Cyy,J(k) X + = 
<fyy, x \(k) X - and Gzz, x \(k) X + = Ozz,x\{k)x- requires second order terms (x - x k ) 2 , while higher-order gradients require the 
higher order contributions omitted in previous equations. The continuity functions appearing in (13), (13)' are 
computed in a straightforward way by enforcing previous conditions. 



III. ENERGY UPDATING TECHNIQUE 

Equivalent forms of models by the energy standpoint can be constructed using energy-based weak form 
versions of governing equations. As shown in [36] and [37], an iterative post-processing technique working on a 
spline interpolation of FEA result can be used to construct an updated solution that locally improves the 
accuracy of a finite element analysis carried out using standard shear-deformable plate elements. As shown in 
[38], a modified expression of the displacements field by a structural model that is free from derivatives of the 
functional d.o.f. can be constructed in order to obtain a CO equivalent model without any additional d.o.f. using 
symbolic calculus to obtain once for all the relations in closed-form. The comparison with exact solutions of 
sample test cases demonstrated that the equivalent model (EM) can provide equally accurate results than the 
original model (OM), even when the through -the -thickness stress distributions are rather intricate. This opens 
the possibility of developing an efficient finite element free from nodal derivatives of the d.o.f. 

Hereafter, the technique developed in [38] is retaken and revised with the aim of constructing a version 
of SEUPT that produces an equivalent zig-zag model suited to develop an efficient and accurate mixed plate 
element for analysis of laminated and sandwich composites. 

3.1 Equivalent zig-zag model 

The basic assumption of SEUPT is that each displacement of the OM model appearing in Eq. (1), 
hereon indicated as u(x, y, z) 0M , v(x, y, z) 0M , w(x, y, z) OM (superscript OM will define all the quantities that refer 
to this model) can be reconstructed using energy equivalence concepts, so that all the derivatives of the 
functional d.o.f. can be replaced. The objective is to derive modified expressions of the displacements, i.e. u(x, 
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u(x, y, z) BM , v(x, y, z) BM , w(x, y, z) BM (superscript EM will refer to the equivalent model) that will be used within 
each energy integral, so to have the energy functional free from derivatives. 

OM 

It is assumed that the displacements of the OM model (in compact form v ) can be recast as the 
sum of a part v free from derivatives of the d.o.f. and terms v containing the derivatives 

\/° M (x, y, z) = \/ 0 (x, y, z) + (x, y, z) 

The basic assumption of SEUPT is that each term V u can be replaced incorporating corrective terms 
that are free from derivatives A u 0 , A v ° , Aw 1 , Ay x , Ay ° (in compact form A V 0 ) 

V EM (x, y, z) =V 0 (x, y,z) + AV 0 (x, y, z) (14') 
whose expressions are derived from the energy balance 

<yj(.) I £ = <?]"(•) \ At -<?}(.) \ A f +Sf(.) I,,„= 0 (15) 

in order to be equivalent by the energy standpoint. The principle of virtual works accounting for inertial 
forces is here used as the energy balance, but any other equivalent canonical functional could be successfully 

employed. As customarily, the symbols [a } and {e .. } represent the strain and stress vectors, respectively, b. 
and t. are the body and surface forces, p is the density and w.are the displacements. The expressions of 
corrective terms A V 0 are specific for any sub-integral constituting the virtual variations of the strain energy 

1 T 

J (.) \ A = —J I cr y j {f ..} dV , work of external forces J (.) \ Af = J b.u.dV + J t.u.dS and work of inertial 

^ V v s 

forces | (.) \ Am = J - p u i a . dV Therefore, each corrective term is computed by splitting the energy balance 

V 

into five independent balance equations, here indicated in compact form as 

[*/<.)!,]" =0; [*/(•> I, ]' =0; [*/(.)!,]" =0; [s J (.) \ E f = 0 ; [ 8 J (.) I E = 0 ; (16) 
i.e. one for each primary variable 8u 0 , 8v° , Sw° , Sy , 8y , and then further splitting into all 

sub-integrals <5j"(.) I £ that added together form the five separate contributions (16) (i.e. 

s<y|(.) I £ _^ = 5\Q \ E ) 

Different corrective terms Am", Av° , Ate", A/ , "are computed for each sub-integral 

8 1 (.) I by imposing the result by the OM model to be equal by the energy standpoint to its counterpart 

by the EM model, in symbols: 



['I 



\y OM 

C) I 



r" 



r ~\r°OM r 



0 1 ] 

£ > J 


u EM 

= 0 ; 




'° EM 


Ol ] 

E > J 


= 0 ; 






oi ■ 


= 0; 


E > - 






y ° EM 


'O I B _^ 


J" =o; 










(•) 1 

E > 


J =o; 



(17) 



The displacement fields are described as Hermite polynomials within the sub-integrals with the 
superscript 0M , while Lagrange polynomials are used within sub- integrals with the superscript EM , because the 
functional d.o.f. and their derivatives should be assumed as nodal d.o.f. whether a conforming element is 
derived by the OM model, while no derivatives are required for an element derived by the EM model, as in the 
current case. At least third order derivatives are involved by the OM model as a consequence of the enforcement 
of the interfacial stress continuity and of the stress boundary conditions, but higher order derivatives could be 
involved enforcing additional conditions. 
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Representing each functional d.o.f. (i.e. u°(x,y), v°(x,y), w°(x,y), y (x,y), y°(x,y))in 
compact form with the symbol p , the inner representation of variables of the OM model is expressed as the 
sum of tensor products of Hermite's polynomials P H_ in x and y of the following form (index i represents the 
nodes numbering): 

4 4 4 4 4 4 

1=1 1=1 i=l :=1 i=l 1=1 

(18) 

Explicit terms have been reported in Eq. (18) only up to the second order of derivation, but higher- 
order terms could be included since the order of the interpolating functions is a direct consequence of the order 
of derivation, which in turns depends on the type of physical constraints enforced. 

1-lEM 
5 1 (.) I J , because the EM model is a C° 

model. Hence, a computationally efficient Lagrangian representation can be used for the EM model, the same 
discussed forward to develop the present finite element. 

The corrective displacements are obtained from Eq. (17) as functions of the nodal unknowns using 
symbolic calculus. Their expressions hold irrespectively of the loading and boundary conditions, lay-up and 
geometry. Being obtained once for all in closed-form, they consistently speed up computations. 

Former balance equations (17) state that since the consistent displacement field by the OM model 
satisfies the energy balance, the modified displacements by the EM model satisfy it too, thus they represent an 
admissible solution by the energy standpoint. Because the EM model is just equivalent form the energy 
standpoint to the OM model, not in a point form, it just provides a correct solution in terms of displacements, 
but it could be unable to satisfy all the boundary constraints in a point-wise sense. In particular, the continuity of 
displacement derivatives and stresses at nodes and sides could not be identically satisfied by the Lagrangian 
representation used within the EM model. On the contrary, it can be satisfied by the Hermitian representation 
used for the OM model and by a suited definition of contributions A L , A c -' p to displacements. Hence, the EM 
model just provides the solution in terms of displacement d.o.f. at any point, while all the quantities that contain 
derivatives, like the stresses, should be determined by the OM model using the analytic expressions obtained 
through symbolic calculus. These operations are indicated in compact form as: 



[<?_[(.) ij" EM \s\Q I J 

[*{(•) 1 J 
[*/<•) I J" 

[•vy EM r -I 



(.) 1 


-in EM 
E J > 




Tv°EM 


O 1 


A -> 


(.) 1 


•m n EM 
E J > 




ly I EM 


(.) 1 


E J ~~ 




y° EM 



; (19) 

y x "OM 
y v °OM 



These energy balance equations could be used as a post-processing technique for obtaining a 
progressively refined solution from the results of a preliminary finite element analysis: 

[ih° EM r -ii/'VjM 

[■|v 0 £'Af r tv OM 

<?{(•) U ^ kjo U ; 

[tw°EM r -iw [, OM 

J/0 ij « Q ij ; (20) 

[ly , EM r -ly , OM 

[■yy ° EM r -iy ° OM 

like in the former applications of SEUPT, but this is not of interest in this paper. Instead, the objective 
is the obtainment of accurate results without any post-processing operation. Indeed, Eqs. (19) and (20) are fast 
carried out using expressions obtained through symbolic calculus, but they however result in an increased 
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computational burden. A mixed element is developed to overcome the use of Eqs. (19), (20) as the stresses, 
which are interpolated separately from displacements, are predicted very accurately. The interlaminar stresses 
that are the ones critical by the viewpoint of strength, stiffness and durability are assumed as nodal d.o.f. 

In order to enable the development of such a mixed element, the interlaminar stresses computed by the 

equivalent displacements of the EM model, here indicated as a J xz , ct J yz a 3 yZi are modified incorporating 

unknown continuity functions A*" 

s 

a' a ' • £ A // (21) 

k = \ 
S 

a J v _ = a J n + £ @ ik) H k (22) 

k =1 
S 

a ' a ' • £ I //. (23) 

* = i 

That make continuous the stresses a J xz , a J yz , <r J zz at the layer interfaces 

= a a v . = ct 2 ct = ct (24) 

The new continuity functions appearing in (21)-(23) are computed in a straightforward way since now 
they are not expressed in terms of displacements as before (9), thus no derivatives are involved, differently from 
[38] where the stresses were obtained from the displacements. Accordingly, their in-plane variation can be 
assumed in Lagrangian form like for the displacement d.o.f. In the present element, the mid-plane stresses are 
assumed as stress d.o.f. Once the finite element solution has been found, their variation across the thickness is 
computed by the constitutive equations (21)-(23) at any z. 

It could be noticed that for keeping the model in C° form, the continuity of the transverse normal stress 
<y 3 Z z,z + >= cr J zz , z <_) appearing in Eq. (11) has been disregarded in Eq. (24). In order to indirectly fulfil this 
necessary condition, the interlaminar stresses of the EM model are made consistent with those of the OM model 
that fulfil (1 1) by further enforcing the equivalence conditions 

r -lOM r -lEM 

[j" (ct B Se K )\ R > J = [5 $ ((ct r _ + A ct t; )Ss R ) l £ ^ J ; 

r -lOM r -\EM 

[J (a w Ss K ) \ £ > J =[S j ((a r + Aa v: )Ss v: )\ e _ j ; (25) 

r -lOM r -lEM 

[j (<J zz Ss a ) l £ _ J = [S j ((a ;: + Aa zz )Ss :t )\ £ > J 

from which the corrective stresses A a , Act , Act are computed in closed-form once for all 

again using symbolic calculus. The Hermitian representation for the OM model and a Lagrangian one for the 
EM model are still used. 



3.2 Finite element 

Based on the updated version of SEUPT discussed above, an eight-node mixed plate element with 
standard Lagrangian interpolating functions is developed. The terms mixed is used to indicate that the master 

fields are internal fields. As nodal d.o.f. the nodal components of displacements u(x, y,z) , v(x, y,z) , 
w(x, y , z) EM and the stresses ct 1 , a \_ , a [_ at the reference middle plane are assumed (see Figure la), as 

mentioned above. In this way, only the continuity of interlaminar stresses is considered. It should be noticed 
that the continuity of stress components parallel to an interface should not be forced, being physically incorrect 
for multi-layered structures. 

Since stresses and displacements can be varied separately, the Hellinger-Reissner multi-field principle: 

Sn „ = 5 [ a s". - —a..S.. u a u dV -off b.u.dV + f t.u.dS ) = 81! „ - SW „„ (26) 

is the governing functional. The symbol Spi represent the components of the elastic compliance tensor, 
the inverse of the elastic stiffness tensor Cyu, so that a ki= £° represents the strains obtained from the 
stress-strain relations 

Accordingly, Vi a y S\^\ a m represents the complementary energy density in term of the master stress 
field. The slave fields are the strains e " and s * : 
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s" „ - C „,,a ,,; s" .. = —(u . + u . .) (27) 

The vector of the nodal d.o.f. for the generic element (e) in its natural reference system 1 is thus 
represented by: 

V} ,\v, EM ) ,{ Wi EM \ Act J] AaJ) A* J] ) (28) 



where , is the node number (z'=l, 8). 

Nevertheless the transverse normal stress gradient is not assumed as a functional d.o.f, in order to 
preserve from stress derivatives as nodal d.o.f, the numerical results will show that the boundary conditions (7) 
will be spontaneously satisfied since many physical constraints can be point-wise enforced, thus it is sufficient 
to enforce remaining ones in a weak form to obtain results in a very good agreement with exact solutions, as 
shown by Icardi and Atzori [48] and Icardi [49] using solid and a layerwise plate elements, respectively. 
Accurate solutions were obtained either for regular problem involving laminates and sandwiches (even when 
they are thick and the properties of their constituent layers are distinctly different), or singular problems like the 
case of a two material wedge with certain side angles and materials with dissimilar properties. Even no case was 
found where the lack of the enforcement of the continuity of the transverse normal stress gradient gave incorrect 
results, it could be noticed that this continuity can be accounted for by a suited definition of contributions (13), 
(13)'. 

Hereon the basic steps towards derivation of the present element will be summarized. Details 
concerning standard aspects will be omitted, being explained in common textbooks. The displacements, 
represented in synthetic form as D and the interlaminar stresses as S, are interpolated separately as: 

8 8 

D = y D e N" r ; S = y S'N " (29) 

/ i disp / i stress * y 

1 1 

since the basic assumption is that the intra-element equilibrium conditions are met in an approximate 
integral form, thus differently to customary mixed elements the nodal stresses are not required to be represented 
in a point-wise self-equilibrating form. Accordingly, the same inner representation in x, y can be assumed both 

for displacements and stresses, i.e. N ' = N " sims =N. In the present case, computationally efficient parabolic 

Lagrange's polynomials are employed. This directly follows form the Hellinger-Reissner governing functional 
(26), since no second or higher-order derivatives of the primary variables (28) are involved, thus no derivatives 
of the d.o.f should be chosen as nodal quantities, consequently C 1 or higher-order interpolating functions are 
unnecessary. 

The option of using the same representation was chosen because it makes easier the development of 
mixed elements without compromising accuracy, as shown since the pioneering papers by Loubignac et al. [39] 
and Nakazawa [40] and since an identical interpolation for displacements and stresses is acceptable by the 
viewpoint of numerical stability of mixed elements, as discussed forward. 

The present element is not intended for use with nearly incompressible materials, because the mean 
stresses, i.e. the pressure, are not treated as independent from the total stresses. However, numerical tests 
demonstrated that accurate results can be obtained in these cases, as shown for a two-material wedge with one 
half incompressible [48]. 

The following standard interpolating functions are used for corner nodes (1, 2, 3, 4) 

= -(i-<r 0 .)(i + c-iry.x^ . •(-lry. -i) (30) 

4 

and at mid-side nodes (5, 6, 7, 8) 

n 5 = -a-'7 oS )(i-£„ 5 2 ) 



^6 = -V-vJ)H+$ 06 ) 



2 

n s = -d-njHi-^A 

2 



(31) 
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This parabolic representation is chosen just in order to obtain accurate results with a relatively coarse 

meshing. 

Representations (30) and (31) meets the compatibility condition, as the displacements and the 
interlaminar stresses are continuous inside and at the edges of the elements. Consistency is also ensured at the 
interfaces. However, since the membrane stresses are not assumed as nodal d.o.f., they are not continuous 
across the interfaces, as mentioned above. But this incompatibility should not be removed since the theory of 
elasticity prescribes that in-plane stresses are discontinuous across the material interfaces. Of course, z should be 
the coordinate in the direction across the layer stack-up, while x, y lie in a plane parallel to the plane of layers. 
With this assumption, the stresses within the finite element are represented as: 

\ \o] } 

(32) 



[ 6 ] being the matrix that defines the in-plane stress components 

K °„ ^/ =[**][*]{<?«,} = p]{<W (33) 

[S *] being the related elastic coefficients and [B ] the derivatives of the shape functions. 

The "stiffness" matrix is derived in a straightforward, standard way from the internal energy functional 



f/ H R substituting the discretized displacements and stresses: 

1 



[ K ] {«(.,} = {«(.,} l v [ 6 ] i B ]-~[ 6 ] l S ][ £ ] dY {*f> 



(34) 



while the vector of generalized nodal forces is obtained from the work of forces W H r substituting the 
finite element discretization of displacements. Assembly of the stiffness matrix and of the vector of nodal 
forces is carried out with the standard techniques. 

The integrals defining the stiffness matrix and the nodal force vector are computed either in exact form 
using symbolic calculus, so to avoid numerical instabilities, or through Gaussian integration, because a reduced 
integration scheme could be used (even if it was never applied for obtaining the numerical results given next). 

To standardize the computation of integrals, as customarily a topological transformation is carried out 
from the physical plane (x, y) to the natural plane {tj, 9) that transforms any quadrilateral element into a square 
element with unit sides 



x.N 



and y = ^ y .N . 



(35) 



In this way, the Jacobian matrix 

r dx dy ] 

| dx ay I 
[dr/ dii J 



[/] 



(36) 



necessary for obtaining the physical derivatives - 



, — from the derivatives 

dx dy 



■ over the 



natural plane 







r d i 


dx 


• = [/]-'< 




d 


d 


, 5 >\ 




d rj 



(37) 



can be computed in an efficient way. As well known, stability, solvability and locking phenomena of 
mixed elements are governed by rather complex mathematical relations (see, Babuska [50] and Brezzi [51]). In 
particular, certain choices of the shape functions could not yield meaningful results. A necessary condition for 
solvability, which in most cases suffices for acceptability of mixed elements, is that the number of displacement 

d.o.f. n must be equal (like in the present element), or larger than the number of stress d.o.f. n . A sufficient 

condition for solvability (see, Olson [52]) is that for a single element that is free of boundary conditions the 
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number of zero eigenvalues of [ K ] e is equal to the number of rigid body modes, the total number of positive 
eigenvalues is equal to the number of stress d.o.f., while the total number of zero and negative eigenvalues is 
equal to the number of displacement d.o.f. However, even when these admissibility tests are passed, erroneous, 
highly oscillating results could be obtained. For homogeneous, isotropic materials this problem is usually 
overcome by relaxing the stress continuity impositions, but in laminated and sandwich composites the 
interfacial stress continuity conditions cannot be skipped, otherwise the physical meaning of solution is lost. The 
opposite possibility is that nevertheless these tests are failed, mixed elements could be effective for solving 
regular problems. Accordingly, mixed elements should be tested considering sample cases with the same 
features of the practical cases to be solved, whose solution in exact or approximate numerical form is available 
for comparisons. 

In the forthcoming section, solvability, convergence and accuracy tests of the element will be 
presented, considering laminates and sandwiches with simply-supported and clamped edges, test cases with 
geometric/material discontinuities and intricate stress variation across the thickness, for which either exact 
solution in closed-form or in numerical form are available in the literature for comparison. 

IV. NUMERICAL APPLICATIONS AND DISCUSSION 

Accuracy, convergence and solvability of the present finite element will be assessed comparing its 
results with exact/numerical solutions of sample test cases available in the literature. Aiming at enhancing the 
intricacies of stress and displacement fields, cases with strongly asymmetrical lay-ups, clamped edges, distinctly 
different properties of layers and high face-to-core stiffness ratios will be analysed. 

For what concerns the through-the-thickness discretization, please notice that all results presented 
throughout this section are obtained considering a number of computation layers equal to the number of physical 
layers, as well as a third order representation of the in-plane displacement and a fourth order representation of 
the transverse displacement for the contributions of Eq. (3). 

4.1 Solvability and convergence 

First of all, an eigenvalue solvability test is presented. In accordance with Olson [52], this test is carried 
out over a single finite element in the shape of a square plate with sides of unit length and free of boundary 
conditions (thickness should be sufficiently small to have a plane-stress, plane-strain problem). The same 
isotropic material analysed by Mijuca [53] in a similar test is here considered with the following mechanical 
properties Ej (= E 2 = E 3 ) =1 and u=0.3. The results of this test are reported in Table 1. 

According to the rules discussed above, the finite element shows a number of zero eigenvalues equal to 
the number of rigid body modes. At the same time, the total number of positive eigenvalues is the same as the 
number of generalized stress d.o.f. and the total number of zero and negative eigenvalues is equal to the number 
of generalized displacement d.o.f. These conditions being satisfied in the test, the sufficient condition for 
solvability is passed by the present mixed element that thus can be applied to solution of other sample cases, in 
order to verify its behaviour in problems of practical interest. 

Aiming at verifying the convergence rate of the present element, the sample case analysed by Zhen et 
al. [23] is considered, with a length-to thickness ratio of 4. The analysis is carried out for a plate that undergoes 
a bisinusoidal loading and is simply supported at the edges. The mechanical properties of the constituent 
material considered for this case are: El/E t =25; Glt/E t =0.5; Gtt/Et=0.2; i) LT =0.2 and the stacking sequence is 
[15°/-15°]. Table 2 reports the through-the-thickness variation of the in-plane stress for different meshing 
schemes by the present element and by Zhen et al. [23], normalised as follows: 

ct I — — , , z I h 

\ 2 2 j 

>' L > (38) 
The results of Table 3 show that the present element is as accurate as the 3D solutions reported in Ref. 
[23], and even with rather coarse meshing the accuracy is satisfactory. This confirms that mixed elements are 
accurate without requiring a very fine discretization as already focused in the literature. 

As a further assessment of the convergence rate, the laminated [0°/90 o /0°] plate for which Vel and 
Batra computed the exact 3D solution in Ref. [42] is considered. The plate is square and has a length to 
thickness ratio (L x /h) of 5. The mechanical properties of the constituent material are the same as in the previous 
case. For what concerns boundary and loading conditions, the plate is simply supported on two opposite edges, 
clamped on the other two and it is undergoing a bi-sinusoidal normal load with intensity p° on the upper face, 
whereas the bottom one is traction free. It could be noticed that the structural model can treat clamped edges 
because a non-vanishing transverse shear resultant can be obtained with all displacement d.o.f. vanishing with 
an appropriate choice of variable kinematics contribution A 1 . The results of Table 3 are normalized according to 
[42] as: 
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10 • H a 



{ 8 ' 2 ' Z J - 100£ 7 -A 3 



w(0,0,z) 

p'l, p i: m 

Like in the previous case, the numbers in square brackets represent the computational time required to 
perform the analysis on a laptop computer with a 1800 GHz double-core processor and 2.96 GB RAM. The 
numerical results show that the present finite element provides very accurate results as compared to the exact 3D 
solution. Also in this case it is shown that even with a rather coarse discretization the accuracy is good. 

4.2 Laminates and sandwiches with different boundary conditions 

Hereafter different cases of analysis available in the literature are considered in order to exhaustively 
verify the accuracy of the present mixed element. 

4.2.1 Simply supported sandwich beam 

Aiming at assessing the capability of present element to smoothly represent the through-the-thickness 
stress and displacement distributions even when material properties of constituent layers abruptly change, it is 
considered an extremely thick, simply supported sandwich beam loaded by a sinusoidal transverse loading. 
Using the present element, the cellular structure of core could be discretized into details, but this is not done 
because the results used for comparisons have been determined considering the sandwich beam as a sandwich- 
like, multi-layered, homogenized structure where the core is described as a quite compliant intermediate thick 
layer and faces as thin, stiff layers. Therefore, it represents for the finite element a laminate having distinctly 
different properties of layers. Of course, this approach can be considered whenever it is supposed that local 
buckling phenomena do not occur in the cellular structure, as they cannot be accounted for by the homogenized 
description adopted. This approach will be used in all the following applications to sandwiches, since it was 
used also to obtain the reference solutions. 

Faces are laminates with stacking sequence (MAT 1/2/3/1/3), where the constituent material have the 
following mechanical properties: MAT 1: E,=E 3 =1 GPa, G 13 =0.2 GPa, u 13 =0.25; MAT 2: E,=33 GPa, E 3 =l 
GPa, G 13 =0.8 GPa, u 13 =0.25; MAT 3: E,=25 GPa, E 3 =l GPa, G 13 =0.5 GPa, u 13 =0.25. The core is made of 
MAT4, whose mechanical properties are as follows: MAT 4: E^E^O.05 GPa, G 13 =0.0217 GPa, u 13 =0.15. The 
length to thickness ratio of the beam is 4 and the thickness ratios of the constituent layers are (0.010h/0.025h 
/0.015h/0.020h/0.030h/0.4h) s . Due to the different thickness of constituent layers and to their distinctly different 
material properties, strong 3D effects rise for this case. However, in order to introduce the further complication 
of a strong unsymmetry, cases are presented with E 3 modulus of the core and upper face layers reduced by a 10 2 
factor. As shown by the exact solution computed in [43] using the technique described in [41], this produces an 
opposite variation of the transverse shear stress across the upper and lower faces that is problematic to capture in 
the simulations. 

The results for this case reported in Figure 2 are normalized as follows: 



(40) 

p p hp hp 

Figure 2 a represents the through-the-thickness variation of the shear stress a „ while Figure 2b 

represents a - , and Figure 2c and 2d represent, respectively, the through-the-thickness variation of the in-plane 

displacement u , and of the transverse displacement w . All the results are obtained considering 150 elements 
over the beam, as shown by the inset reported in Figure 2. In this case, the analysis requires 16 s. The results of 
Figure 2 show that the present mixed element can obtain results as accurate as the exact 3D solutions even when 
the constituent material have abruptly changing mechanical properties that determine strong asymmetrical stress 
and displacement fields, using a reasonably fine mesh. 

4.2.2 Simply supported sandwich plate 

Still aiming at verifying the accuracy of the present mixed element in describing stress and 
displacement fields of structure with strong asymmetric effects, it is now considered the sandwich plate for 
which Brischetto et al. [44] compute the exact 3D solution. In this case, the asymmetry is due to differences 
between the thickness ratios of upper (symbol us) and lower (symbol Is) faces, which are respectively hi s =h/10; 
h us =2h/10 with respect to the thickness h of the plate. The thickness ratio of the core (symbol c) is h c =7h/10. 
The constituent material has the following mechanical properties: Ei s /E lls =5/4, E] S /E C =10 5 , Vi s = v us =v c =v= 0.34. 
For what concerns loading and boundary conditions, the plate is simply supported and it is undergoing bi- 
sinusoidal loading. The plate has a length to thickness ratio L x /h=4 and a length-side ratio L y /L x =3. Figure 3 
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reports the comparison between the exact solutions presented in Ref. [44] and the numerical results for what 
concerns the in-plane displacement and the transverse shear stress computed using the present finite element, 
normalized as follows: 

( L } } ( L \ 

a | 0,— ,z [h u 0,— ,z LB h 

~ - I 2 } " I 2 >— (An 
a „ = ; u = — (41) 

The results of Figure 3 have been obtained considering 300 elements and modelling only a quarter of 
the plate (see inset in Figure 3). Also in this case, despite the asymmetry of the structure, the present mixed 
element can obtain results in a very good agreement with the exact 3D solution at any point with the right 
gradients at the interfaces, requiring just 35 s to perform the analysis. 

4.2.3 Cantilever sandwich beam 

In order to more deeply assess whether the present finite element can accurately describe the stress 
field also with clamped edges, the cantilevered sandwich beam subjected to a uniform transverse loading of Ref. 
[45] is considered. The beam has laminated faces made of unidirectional Carbon-Epoxy laminates (E,= 157.9 
GPa, E 2 =E 3 = 9.584 GPa, G 12 = G 13 = 5.930 GPa, G 23 = 3.277 GPa, v 12 = v 13 = 0.32, v 23 = 0.49), while the core is 
made with a PVC foam core (E= 0.1040 GPa, v= 0.3). According to [45], a length-to-thickness ratio of 10 is 
assumed, while the thickness ratios of the constituent layers are 0. lh/0.8h/0. lh. 

Figure 4 reports the through-the-thickness distribution of the transverse shear stress and the meshing 
scheme adopted. The stress is normalised as follows: 

— 2h ( L L "| 
L x p° \ 5 2 ) 

The results of Figure 4 confirm the accuracy of the present mixed element as it obtains results in good 
agreement with those by the 3D finite element of Ref. [45]. For what concerns the computational times, only 
17.5 s are required to perform the analysis. Thus, even if the model of Eqs. (1) and its EM counterpart have as 
functional d.o.f. the displacements and the transverse shear rotations at the middle plane, it does not give poor 
results for clamped structures like the equivalent single-layer models with classical five mid-plane d.o.f., 
because enforcing the vanishing of displacements at the clamped edges does not necessarily mean that the 
transverse shear stress resultant also vanishes. Also in this case, quite accurate results are obtained with a not 
extremely refined meshing. 

4.2.4 Cantilever piezo-actuated beam 

With the goal of verifying the accuracy of the present mixed element for clamped structure with 
abruptly changing material properties, the cantilever piezoactuated beam formerly studied by Robbins and 
Reddy [46] is considered. The structure is made of an underlying aluminium beam substructure (thickness 15.2 
mm), an adhesive film (thickness 0.254 mm) and a piezoactuator (thickness 1.52) bonded on the upper face. The 
material properties of the constituent materials are: E a i um =69 GPa, G a i llm =27.5 GPa, u a i um =0.25; E ad =6.9 GPa, 
G ad =2.5 GPa, u ad =0.4; E lpiezo =69 GPa, E 3piezo =48 GPa, G piezo =21 GPa, T)i 3p ie Z o=0.25, u 31piezo =0.175. Piezoactuator, 
adhesive and aluminium substructure form an unsymmetrical three material laminate that exhibits 
bending/extension coupling. The only acting loads for this case are the self -equilibrating loads induced by the 
piezoactuator. A bending deformation is provided by applying an actuation strain of 0.001 to the piezoelectric 
layer via an applied electric field. The length of the beam is 152 mm. 

Figure 5a shows the spanwise variation of membrane, transverse shear and transverse normal stresses 
near the top of the aluminium substrate (zl= 6.61 mm) and at the centre of the adhesive layer (z2= 6.84 mm), 
while Figure 5b represents the stress fields across the thickness close to the free edge. The stress distributions 
are presented in the following normalized form, according to [46], where the exact solution for this case is 
presented: 

- A. -a -10 3 
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the subscripts alum, ad, piezo and tot being used to indicate the cross sectional area and the elastic 
moduli of the substrate structure, of the adhesive and of the piezoactuator layer, respectively, which are 
considered as isotropic materials. 

The results of Figure 5 show that the present mixed element is accurate also for this case in which a 
singular-like stress filed arise. In details, in Figure 5a it can be noticed that the element accurately describes the 
unwanted dangerous stress concentrations that appear at the free end of the beam that could cause debonding of 
the piezoactuator in service, without showing erroneous oscillating results close to the edge. Also the through - 
the-thickness distribution of the stresses is well described as shown in Figure 5b, in particular the mixed element 
is capable to capture the sharp variation of the bending stress that can be seen near the bonding layer also in this 
case with a reasonably refined meshing. For this case of analysis the present element required 47 s. 

4.2.5 Double lap joint with aluminum adherents 

As an application to a practical case, it is now considered the joint with aluminium 2024 T3 adherents 
analysed by Nemes and Lachaud [47]. The adhesive is epoxy resin REDUX 312/5 (E=27 GPa, G=l GPa, 
t)=0.35) and it is 0.1 mm thick. The outer layer is 2 mm thick, the inner one is 4 mm thick, and the overlap 
length is 50 mm. A pressure load in x with intensity 1 N/mm is applied to the inner adherent, while the lower 
and the upper adherents are clamped at right and left bounds. 

This case is treated as a plate with step variable properties in the finite element analysis, the adherents 
being constituted by a homogeneous and isotropic medium, while the overlap is constituted by and a three-layer 
laminate. Of course, all continuity functions vanish in the adherents, while they exist at the two interface 
between adhesive and lower adherent, adhesive and upper adherent. The piecewise contributions of Eqs. (13) - 
(13') are considered in this case in order to make continuous the stresses and their gradients at the interface of 
adherents and overlap, where properties vary from a single layer to a three-layer plate. Specifically, the terms of 
first order are aimed at fulfilling the continuity of the stresses; those of second order are aimed at satisfying the 
continuity of the stress gradients. The variable kinematic contributions across the thickness (3) in this case are 

aimed at satisfying also the stress boundary conditions at the ends of the overlap Q p CI : 

N = f a dz = - f a _dx; M = f z a dz; Q = f a dz; (44) 

Jn, " Jn, " Jn, Jn, x * 

N, M, Q being the in-plane, bending and shear resultants, respectively, along with the stress -free 
boundary conditions 

a = 0; f a dz 
a =0; f a dz 

XX J Q " 

at the free edges of the overlap. The results for this case are presented in Figure 6, where the in-plane 
variations of the peel stress and of the shear stress are reported. They are obtained requiring 69 s to perform the 
analysis and modelling the whole joint with 600 elements with 2 concentrations of elements across the variation 
of adherents, as shown by the inset. Also this case confirms the accuracy of the mixed element, which correctly 
describes the stress peaks near the end of the overlap, which can have detrimental effects on the service life of 
the joint. 

V. CONCLUDING REMARKS 

A quadrilateral, eight-node, mixed plate element for analysis of laminates with general lay-up and 
sandwich composites with high face-to-core stiffness ratios and distinctly different properties of top and bottom 
faces, was developed. It has the displacements and the interlaminar stresses at the reference mid-plane as nodal 
d.o.f. Its structural model is a 3D physically based zig-zag one with fixed functional d.o.f. (the three 
displacements and the two transverse shear rotations at the mid-plane), whose representation can be refined 
across the thickness without increasing the unknowns. This model also allows the material properties to 
suddenly vary in the in-plane directions. A piecewise variable transverse displacement is assumed, since it has a 
significant bearing for certain problems, e.g. when temperature gradients cause thermal stresses, for the crushing 
behaviour of sandwiches, around geometric and material discontinuities and close to loading points. The 
expressions of continuity functions and higher-order coefficients of displacements incorporated in order to a 
priori satisfy all the physical constraints are determined once for all in closed-form as expressions of the five 
d.o.f. and of their spatial derivatives using a symbolic calculus tool. 

The energy updating technique (SEUPT) is used to convert the d.o.f. derivatives without resulting in an 
increased number of variables. It allows obtainment of an equivalent C° version of the zig-zag model that is used 
to develop the finite element. The three displacements and the three interlaminar stresses at the mid-plane are 
assumed as nodal d.o.f. All these quantities are interpolated using the same standard serendipity shape functions 



= 0; £ J (X xz dz = 0 ' 

^ k (45) 
= 0; £ J a tt dz = 0; 
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because this simplifies the development of the element without resulting in any accuracy loss, as shown by the 
comparison with exact solutions of sample cases. Consequently, equilibrium equations are satisfied just in 
approximate integral form. 

Accuracy, solvability and convergence were assessed considering multi-layered monolithic structures 
with cross-ply and angle-ply lay-up and sandwich-like structures with also abruptly changing properties of 
faces, with different loading and boundary conditions, for which exact solutions are available. The structural 
model and the mixed formulation make easier the enforcement of stress -boundary conditions, either because the 
coefficients of higher-order terms can be determined from enforcement of these conditions, or the stresses are 
nodal variables. An application was presented also to an adhesively bonded joint, which represents a case with 
step varying properties, the adherents and the overlap being seen as interfaced plates with a different lay-up. The 
element was shown to be fast convergent and its results were not affected by the mesh configuration, as accurate 
predictions were obtained also with rather coarse meshes. It appears as a good alternative to the ones to date 
available in the literature, as it offers a good accuracy with the minimal number of unknown nodal variables at 
an affordable cost. 
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Table 1. Solvability test for an isotropic material 
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Table 2. Normalised in plane stress for a [15°/-15°] square plate by Zhen et al. [23]and by the present mixed 

element with progressively refining meshing. 



www.ijceronline.com 



Open Access Journal 



Page 21 



(f mixed layerwise quadrilateral plate element with variable out-of- plane kinematics and fixed d.o.f. 



z/h 


-0.50 


-0.40 


-0.30 


-0.20 


-0.10 


0.00 


0.10 


0.20 


0.30 


0.40 


0.50 




Ref. [42] 


0.000 


2.609 


3.227 


2.643 


2.113 


2.093 


2.100 


2.668 


3.340 


2.734 


0.000 




Present (3x3) [1.97] 


0.000 


2.720 


3.830 


2.949 


2.450 


2.360 


2.450 


2.950 


3.950 


2.920 


0.000 


a u 


Present (5x5) [2.70] 


0.000 


2.650 


3.450 


2.734 


2.325 


2.248 


2.299 


2.750 


3.430 


2.805 


0.000 




Present (9x9) [9.26] 


0.000 


2.607 


3.227 


2.643 


2.114 


2.094 


2.099 


2.668 


3.340 


2.735 


0.000 




Ref. [42] 


1.152 


1.158 


1.162 


1.167 


1.173 


1.180 


1.188 


1.198 


1.208 


1.218 


1.227 




Present (3x3) 


1.613 


1.545 


1.450 


1.406 


1.375 


1.340 


1.333 


1.388 


1.398 


1.415 


1.443 


w 


Present (5x5) 


1.234 


1.288 


1.291 


1.299 


1.304 


1.306 


1.309 


1.295 


1.299 


1.305 


1.312 




Present (9x9) 


1.153 


1.158 


1.161 


1.168 


1.172 


1.180 


1.188 


1.197 


1.207 


1.217 


1.228 



Table 3. Normalised shear stress and transverse displacement for a [0°/90 o /0°] square plate by Vel and Batra 



[42] (exact 3D solution) and by the present mixed element with progressively refining meshing. 
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Figure 1. Reference system for the element in: a) the natural plane and b) the physical plane. 
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(f mixed layerwise quadrilateral plate element with variable out-of- plane kinematics and fixed d.o.f. 
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Figure 2. Stress and displacement field by the present mixed element and exat 3D solution [43] for a sandwich 

beam. 
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Figure 3. Transverse shear stress and in-plane displacement by the present mixed element and by Brischetto et 
al. [44] (Exact solution) for a simply supported sandwich plate. 
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(f mixed layerwise quadrilateral plate element with variable out-of- plane kinematics and fixed d.o.f. 
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(f mixed layerwise quadrilateral plate element with variable out-of- plane kinematics and fixed d.o.f. 
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Figure 6. Span-wise distribution of a) shear stress and b) peel stress by Nemes and Lachaud [47] and by the 

present model 
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